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1 .  INTRODUCTION 


In  1981-1983  a  computer  program  was  developed  which  solves  numerically 
two-dimensional  Euler  equations  using  the  Godonov  method  (Ref.  1).  The 
program  was  intended  primarily  for  turbomachinery  application,  but  with 
modification  to  the  boundary  conditions,  application  to  any  problem  which  is 
modeled  by  the  Euler  equations  is  possible. 

This  report  covers  one  year  of  extensive  testing  of  the  code  for 
different  flow  regimes.  This  testing  was  done  with  three  main  objectives: 

a)  To  gain  experience  in  simulation  of  the  various  problems  related  to 
the  research  activity  in  the  Turbopropulsion  Laboratory, 

Naval  Postgraduate  School,  Monterey,  California. 

b)  To  test  code  accuracy  by  comparing  numerical  with  analytical 
solutions. 

c)  To  find  the  source  of  errors  for  various  numerical  simulations  and  the 
basic  Godunov  method  in  order  to  reduce  these  errors. 

Some  of  the  results  of  these  experiments  with  the  code  were  reported  in 
other  papers  and  appear  as  appendices  to  this  report.  Other  cases  were 
simulated,  but  not  published,  and  will  be  reported  here.  Some  unsuccessful 
attempts  to  improve  the  code  accuracy  are  also  described. 

The  topics  studied  and  the  progress  made  are  described  in  the  following 
sections. 

1)  The  problem  of  the  gradual  opening  of  wave  rotor  passage. 

2)  Numerical  modeling  of  the  nonsteady  thrust  produced  by  intermittent 
pressure  rise  in  a  diverging  channel. 

3)  Numerical  techniques  for  wave  rotor  cycle  analysis  (A.  Mathur). 

4)  Development  of  a  grid  generation  program  in  order  to  generate  a  grid 
over  wedge-arc  configuration. 

Test  case  of  the  supersonic  flow  in  a  channel  with  arc,  wedge, 
wedge-arc,  sinsoidal  arc-bumps, 

6)  Modification  of  the  Godunov  method  to  include  shock  fitting  using  the 
oblique  shock  wave  theory. 

7)  Development  of  the  one-dimensional  code  based  on  Verhoff's  method. 

8)  Modification  of  the  program  for  cylindrical  geometry  problems,  and 
modeling  of  the  flow  in  CDTD. 

9)  Modeling  of  the  flow  in  a  detonation  engine, 

10)  Development  of  the  SELCO  Method. 

11)  Grid  generation  for  CDTD  inlet 


PROBLEMS  CONSIDERED 


2.1  The  Gradual  Opening  in  Wave  Rotor  Passage. 


Influence  on  the  flow  pattern  of  the  gradual  passage  opening  in  the  wave 
rotor  was  studied  on  the  basis  of  a  numerical  simulation.  It  was  found  that, 
in  most  cases,  a  significant  volume  of  the  passage  will  have  rotational  flow, 
which  should  lead  to  the  mixing  between  the  driver  and  driven  gases.  In  some 
cases,  losses  will  occur  as  a  result  of  the  multiple  reflections  of  the  shock 
and  pressure  waves  from  the  passage  walls.  It  was  shown  that  the  interface 
between  driven  and  driver  gas  will  be  oblique  to  the  passage  walls  when  the 
passage  opens  gradually,  and  that  the  Interface  can  retain  its  obliqueness  to 
the  walls.  The  results  for  the  rectangular  and  skewed  passages  are  reported 
in  Appendix  A. 

In  the  case  of  the  skewed  passage,  the  perfect  compression  of  the  driven 
gas  with  very  small  mixing  losses  could  be  achieved  if  the  channel  opens 
according  to  a  certain  formula  of  the  optimum  opening  velocity.  It  was  found 
that  the  velocity  of  opening  could  be  predicted  by  the  equation: 

V  ^sh  , , V 

''open  =  j  (1) 

Sin  uqi* 


where:  ''open  =  Che  velocity  of  the  gradual  opening 

'^sh  =  the  shock  wave  velocity  in  the  driven  gas 

Lsk  “  angle  of  the  skewed  passage 


2.2  Numerical  Modeling  of  the  Nonsteady  Thrust  Produced  by  Intermittent 


Pressure  Rise  in  a  Diverging  Cnanne 


The  dynamics  of  the  expansion  of  detonation  products  in  a  diverging 
channel  were  Investigated  numerically.  The  influence  of  particular  features 
in  the  expansion  process,  such  as  the  presence  of  reversed  flow  and 
propagation  of  hammer  shocks  on  the  production  of  thrust  were  examined. 
Sequential  expansions  of  detonation  products  were  also  modeled  and  it  was 
concluded  that  in  order  to  maintain  a  high  frequency  periodic  mode  of 
operation  for  propulsion  applications,  the  channel  should  be  refilled  with 
ambient  air  after  each  expansion.  The  influence  of  the  ratio  of  ambient  air 
to  detonation  products  volume  on  the  dynamics  of  the  thrust  production  and  on 
the  impulse  generated  during  the  expansion  are  also  reported.  The 
presentation  of  the  results  of  this  investigation  is  given  in  Appendix  B. 

2.3  Numerical  Techniques  for  Wave  Rotor  Cycle  Analysis. 


With  the  author's  participation,  A.  Mathur  developed  a  one-dimensional 
code  based  on  Godonov's  method  for  wave  rotor  cycle  analysis.  Because  the 
resolution  of  the  code  on  the  interfaces  was  poor,  it  was  decided  to  develop  a 
one-dimensional  code  based  on  the  random  choice  method.  Non-standard  boundarv 
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conditions  for  this  problem  were  developed.  A  paper  was  written  and  accepted 
for  publication  in  the  proceedings  of  the  ASME  Symposium  on  Nonsteady  Flows. 
The  paper  is  attached  as  Appendix  C, 


The  grid  generation  routine,  which  was  used  for  the  cascade  and  wave 
rotor  simulations  did  not  allow  a  full  control  of  the  grid  parameters.  This 
control  was  especially  needed  for  the  simulation  of  the  test  problems,  since 
we  wanted  to  insure  that  the  grids  have  the  same  general  parameters  for  the 
different  wedge-arc  configurations  considered. 

A  program  which  allows  the  grid  generation  with  full  control  over  grid 
parameters  was  developed.  The  program,  with  an  example  of  Che  grid,  is  given 
in  Appendix  D  was  delivered  to  the  Naval  Postgraduate  School  and  stored  in  an 
archive  file. 

2.5  Test  Case  of  the  Supersonic  Flow  in  a  Channel  with  Arc,  Wedge,  Wedge-Arc, 
Sinusoidal  Arc-Bumps. 


Performance  of  the  Godunov  code  was  tested  for  numerous  supersonic  and 
subsonic  Internal  flow  problems.  The  tasks  included: 

a)  An  attempt  to  determine  what  causes  the  nonsymetrical  behavior  for 
subsonic  flow  case  and  whether  it  is  possible  to  increase  the 
accuracy  of  the  basic  Godunov  code  by  changing  boundary  conditions. 

b)  A  comparison  of  the  numerical  solution  with  analytical  predictions 
for  supersonic  flow. 

For  subsonic  flows  it  was  determined  that  the  main  source  of  errors  was  slop 
discontinuity  at  the  corners  of  the  "bump".  When  the  circular  arc  bump  was 
replaced  by  a  sinusoidal  arc  "bump",  the  symmetry  of  the  subsonic  solution 
improved  substantially.  This  resulted  from  the  fact  that  in  the  case  of  the 
sinusoidal  arc  only  the  second  derivatives  are  discontinuous  at  the 
surface.  It  was  found  that  the  boundary  conditions  at  the  surface,  which  are 
used  by  the  basic  Godunov  method  (solution  of  the  Riemann  problem  between  the 
physical  point  and  it  mirror  image  with  the  same  pressure  and  density  but  with 
the  negative  velocity),  produce  the  most  consistent  results  for  various 
problems.  Results  for  this  simulation  are  shown  in  Figure  1.  The  grid  is 
shown  in  Figure  2. 

Efforts  to  Improve  the  subsonic  flow  solution  were  discontinued  when  it 
was  realized  that  a  more  objective  criteria  than  symmetry  of  solution  was 
needed.  For  this  reason,  work  on  supersonic  flows  was  begun,  since,  for  these 
conditions,  some  analytical  solutions  are  available.  First,  the  method  for 
supersonic  flow  (M  »  2.2)  in  a  channel  with  a  10%  chick  circular  arc  "bump  was 
tried.  The  grid  for  this  case  is  shown  in  Figure  3.  In  Figure  4  we  show 
results  for  this  geometry  when  the  inlet  Mach  number  is  2.0.  Results  are 
shown  in  the  form  of  the  distributions  of  the  pressure  coefficient  over  the 
surface  of  the  lower  wall  of  the  channel.  Numerical  results  correspond  very 
well  Co  the  analytical  in  this  case.  It  was  found,  however,  that  when  other 
flow  parameters  are  compared  with  the  analytical  solution,  there  was 
substantial  disagreement  between  the  resits.  Figure  5  shows  the  comparison  of 
the  pressure  coefficient,  entropy,  density,  velocity  and  Mach  number  that  were 


calculated  numerically  with  the  analytical  solution.  The  comparison  is  done 
for  supersonic  flow  with  M*1.6  over  5%  thick  wedge-arc  "bump".  As  shown  in 
Figure  5,  the  errors  are  produced  mostly  on  the  shock  waves  and  that  the 
largest  is  the  relative  error  in  entropy.  The  fact  that  only  pressure  is 
simulated  accurately  was  confirmed  for  the  numerous  test  cases  in  which  the 
arc  thickness  and  the  wedge  angle  varied  in  the  wide  range  of  values. 

2.6  Modification  of  the  Godunov  Method  based  on  the  Shock  Fitting  using  the 
Oblique  Shock  Wave  Theory. 

A  shock  fitting  was  attempted  to  reduce  the  errors  on  the  oblique  shocks 
The  procedure  for  fitting  was  as  follows: 

1)  Calculate  the  values  using  Godunov's  method. 

2)  Using  oblique  shock  wave  theory,  calculate  the  exact  value  of  the 
parameters.  This  calculation  is  based  on  the  assumption  that  the 
pressure  is  computed  correctly  by  the  Godunov  method. 

3)  Substituting  values  behind  the  shock  for  the  exact  values.  It  was 
found  that  the  procedure  did  not  work  if  it  was  applied  as 
described  above.  Many  modifications  of  this  prodedure  were  tried. 
Some  involved  interpolation  of  the  values  ahead  and  behind  the 
shock.  The  last  approach  did  improve  the  accuracy  of  the  entropy 
calculation  but  introduced  oscillations. 


2.7  Development  of  the  One-Dimensional  Code  based  on  the  Verhoff  Method. 

A  new  formulation  of  the  Euler  equations  and  numerical  method  to  solve 
them  were  proposed  by  Verhoff  (Ref.  2).  A  basic  code  implementing  his  method 
for  one-dimensional  flow  was  developed.  In  Figure  6,  results  are  shown  for 
the  case  of  two  colliding  shock  waves  in  a  tube.  These  results  are  a 
reconstruction  of  the  results  presented  in  Ref.  1  in  Figure  5.  Comparison  of 
these  two  figures  shows  that  our  code  implements  the  basic  Verhoff  method 
correctly.  The  listing  of  the  program  is  presented  in  Appendix  E. 

We  were  unable  to  implement  the  shock  fitting  in  our  code.  Comparison 
between  the  analytical  solution  and  numerical  was  10-15%  in  error  for 
velocity,  pressure  and  density  in  the  region  of  the  shock  wave,  if  the  shock 
was  not  fitted. 

2.8  Modification  of  the  Program  for  Cylindrical  Geometry  Problems,  and 

Modeling  ot  the  Flows  in  cDTITi 

The  Godunov  code  was  modified  for  simulation  of  the  flow  with 
cylindrical  velocity.  The  modified  program  was  used  for  modeling  of  the  flow 
field  in  CDTD  device  (Ref.  3).  The  task  of  the  modeling  was  to  test  how  the 
rotational  component  of  velocity  will  influence  flow  field. 

In  Figure  7,  one  of  the  geometries  of  the  CDTD  is  shown  with  a 
computational  grid  covering  the  domain  of  integration.  In  Figure  8,  the 
Iso-Mach  lines  for  the  geometry  are  shown  for  the  inlet  flow  conditions 
indicated  the  figure  7. 


t*  / 


The  conclusions  of  preliminary  research  was  that  the  program  models  the 
flow  correctly  and  that,  for  the  tested  conditions , the  results  are  close  to 
those  given  by  the  program  developed  by  Hirsch  as  discussed  in  Ref,  A. 

2.9  Modeling  of  the  Flow  in  Detonation  Engine. 

The  flow  field  in  the  projected  combustion  chamber  of  a  detonation 
engine  was  simulated  numerically  using  the  Godunov  code.  The  task  of  this 
research  was  to  study  the  dynamics  of  the  discharge  of  the  high  pressure  gas 
after  detonation.  A  number  of  configurations  were  tried.  In  Figure  9,  an 
example  of  this  type  of  simulation  can  be  seen.  Pressure  in  the  detonation 
chamber  drops  from  15  atm  to  0.55  atm  in  0.00131  second.  It  was  important  to 
establish  that  time  since  it  is  one  of  Che  limitations  on  the  detonation 
frequency.  In  this  example,  the  ability  of  the  program  to  simulate  large 
gradients  of  pressure  can  be  seen. 

2.10  Development  of  the  SELCO  Method. 

A  new  method  for  supersonic  flow  simulation  was  developed.  A  full 
description  of  the  method,  with  illustrations  of  its  application  to  the 
supersonic  flow  over  the  wedge,  is  given  in  Appendix  D.  This  method 
demonstrates  that  inaccurate  modeling  of  the  oblique  shock  waves  produced  by 
the  Godunov  method  is  the  result  of  the  obliqueness  of  the  shock  wave  with 
respect  to  the  edges  of  the  cells  of  the  computational  grid  covering  domain  of 
integration.  It  was  shown  that  only  when  the  shock  surface  is  parallel  to  the 
two  opposite  edges  of  the  cell  that  the  oblique  shock  can  be  accurately 
calculated. 
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A  new  method  of  the  Local  Cell  Orientation,  -  SELCO,  was  proposed  to 
allow  local  reorientation  of  Che  cells  in  the  vicinity  of  the  shock  waves. 

The  efficiency  of  the  SELCO  method  was  demonstrated  for  the  simulation  of  the 
oblique  shock  waves  in  the  supersonic  flow  over  the  wedge. 

2.11  Grid  Generation  for  CDTD  Inlet. 


For  a  given  distribution  of  the  pressure  measurement  ports  on  the  CDTD 
blade  surface,  a  computational  grid  was  developed.  Several  variants  of  the 
CDTD  geometry  were  implemented.  However  all  grids  which  were  generated  proved 
to  lead  to  large  computational  errors.  This  work  was  not  continued  due  to  time 
constraints  on  the  project. 

3 .  Conclusions 

The  principal  findings  of  the  program  of  code  testing  were  as  follows: 

1,  The  code  has  a  robust  ability  to  analyse  unsteady  time  dependent 

flows  and  provide  qualitative  information  on  complex  shock  patterns 
and  pressure  fields. 


With  attention  to  boundary  conditions  and  particularly  grid 
structure  it  is  possible  to  reproduce  analytical  flow  solutions  witli 
great  accuracy. 
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3.  The  solutions  produced  are  sensitive  to  grid  orientation.  The  grid 
must  be  re-orientated  locally  along  shock  or  discontinuity  planes. 


4.  Care  must  be  used  in  selecting  the  grid  distribution  on  areas  of 
rapid  curvature  change. 

Despite  the  program  sensitivity  under  certain  conditions,  with  skilled 
use  it  can  be  most  instructive  in  establishing  wave  and  shock  structures 
in  complex  flow  situations.  A  significant  example  is  the  unsteady 
boundary  of  Section  2.1  where  in  it  was  possible  to  generate  a  minimum 
loss  passage  angle  for  a  wave  rotor  device. 
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The  Problem  of  Gradual  Opening  in  Wave  Rolor  Passages 


S.  Eidelman* 

Naval  Postgraduate  School,  Monterey,  California 


1  hr  influrncr  on  Ihr  flow  piurrn  of  Ihr  Krndiiil  pirsagc  oprning  in  thr  wivr  rotor  i«  sludird  on  Iht  basis  of 
numeriral  simtilalion.  It  is  found  that,  in  most  cases,  significanl  volumt  of  Ihr  passage  will  ha>r  rotational  flow, 
which  should  lead  to  the  mising  between  Ihr  driver  and  driven  gases.  In  some  rases,  losses  will  occur  also  as  a 
result  III  the  multiple  reflection  of  Ihr  shock  and  pressure  waves  from  Ihe  passage  walls.  Ii  is  shown  lhal  Ihe 
inlrrfacc  hriwern  driven  and  driver  gas  will  he  <ihliqur  In  Ihe  passage  walls,  when  Ihr  passage  opens  graduallv, 
and  Ihe  interface  can  retain  its  obliqueness  lo  ihe  walls. 


Intrnduclinn 

WAVE  colors  are  devices  that  use  wave  propagation 
ihrough  the  fluid  in  a  rotor  passage  to  transfer  energy 
from  one  fluid  to  another'  or  to  transfer  energy  from  one 
fluid  lo  rotor  shaft  attd  atiother  fluid.’  Wave  rotors,  wave 
engines,  wave  pressure  exchangers,  wave  equalizers,  and 
Comprex'  are  all  based  on  the  same  idea  of  energy  exchange 
in  the  unsteady  waves,  attd  the  topic  of  the  present  study  is 
relevant  to  all  of  these  devices. 

Wave  rotors  offer  the  potential  for  s^nificam  im¬ 
provements  in  air-breathing  engitte  propulsion  cycles  because 
they  can  be  self-coolittg,  siitce  both  high-pressure  gas  and  cold 
low-pressure  gas  use  the  same  passage  for  alternate  periods  of 
lime  in  the  cycle.  Combining  a  wave  rotor  with  conventional 
turbomachincry  components  shows  promise  of  significant 
reduction  in  specific  fuel  consumption  without  weight  or  si/e 
penalties. 

The  principles  of  operation  of  wave  rotor  devices  and  their 
commerical  applications  can  be  found  in  Refs.  1-3.  For 
completeness,  the  general  scheme  of  a  wave  pressure  ex¬ 
changer  will  be  described,  as  illustrated  in  Fig.  1.  One  gas 
(driver)  at  high  pressure  is  used  to  compress  a  second  gas 
(driven).  Ihe  process  is  arranged  to  occur  in  iiibe-like 
passages  wnh  irape/oidal  cross  section  located  oti  the  per¬ 
iphery  of  a  rotating  drum  or  rotor.  The  compression  is 
achieved  sticcessively  in  each  rotor  passage  by  means  of 
compression  waves  or  shock  waves  generated  by  the  entertng 
driver  gas.  1  he  compressed  driven  gas  is  drawn  off  frivm  the 
end  ot  the  passage  when  it  aligns  with  an  outlet  port  The 
driver  gas  iheii  undergoes  a  scries  of  expansions  to  a  lower 
pressure  and  is  scavenged  out  by  freshly  inducted  driven  gas  at 
approximately  the  same  pressure  level.  This  fresh  "charge"  is 
then  compressed  by  the  high-pressure  driver  gas  and  the  cycle 
repeals  iiscll  Steady  rotation  of  the  drum  sequences  the  ends 
of  the  passages  past  stationary  inlet  ports,  outlet  pvirts,  and 
endwalls.  This  establishes  unsteady  but  periodic  flow  pro 
cesses  wntim  the  rotating  passages  and  essentially  steady  flow 
in  the  inlet  and  millet  pons.  The  passage  may  be  oriented 
axially  as  in  l  ig.  I  or  at  a  stagger  angle.  In  general,  wave 
machines  used  as  pure  pressure  exchangers  (e.g.,  "Com¬ 
prex"')  have  usually  axially  oriented  passages,  while  tliose 
with  staggered  passages  may  drive  a  shaft,  sitice  shaft  work 
I'xiraciioii  IS  possible  w  nli  this  latter  confignraiioii. 
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In  the  design  of  the  wave  rotor,  it  is  very  desirable  to 
determine  the  optimal  ratio  between  the  width  and  length  oi 
the  single  passage  of  the  rotor.  Usually  only  two  parameters 
are  evaluated  to  determine  this  ratio:  skin  friciioti  losses  and 
bypass  losses.  The  number  of  passages  on  the  rotor  should  be 
minimal  to  minimize  skin  friction  losses.  On  the  other  hand, 
the  passage  should  be  narrow  compared  to  the  port  widths  to 
reduce  flow  bypass  between  inlet  or  outlet  ports.  The  transient 
process  of  the  passage  opening  or  closing  (as  the  passage  end 
moves  across  a  port  or  moves  from  a  port  to  be  closed  by  the 
endwall,  respectively)  usually  is  not  considered  in  the  design. 
It  is  generally  assumed  that  the  passage  opening  or  cloMire 
occurs  instantaneously . 

Pearson’  tried  to  take  into  account  the  gradual  opening  ol 
the  passage,  assuming  that  the  air  in  the  passage  is  comprev  -ed 
in  a  series  (usually  ilirec)  of  discrete  compression  waves  w  iiidi 
converge  and  ultimately  merge  to  form  a  sluvek  wave.  This 
allowed  him  to  design  a  complicated  wave  machine  cycle  Ut  j 
rotor  using  relatively  short  passages.  However,  since  ihe 
(echnique  was  one-dimensional  it  could  not  reveal  ilie 
peculiarities  of  this  essentially  two-dimensional  How  und 
would  be  valid  only  for  very  weak  waves. 

In  the  present  study,  by  means  of  numerical  modeling,  we 
will  examine  in  real-time  how  the  gradual  opening  mlluei  vcv 
the  wave  formation  in  the  wave  rotor  passages  and  how  ilun 
should  aflect  the  rotor  design. 

Model 

I  hc  assumptions  involved  in  the  numerical  simulations  ,iie 
described  as  follows. 

The  flow  111  each  passage  of  the  wave  rotor  is  unsteadv  and 
periodic.  At  the  same  lime,  the  How  ihroueli  the  pori  o 
(ideally)  steady.'  ’  The  peripheral  width  ol  the  port  is  ummIK 
equal  lo  several  widllis  ot  a  single  passage,  and  herein  n  i- 
avsuiiicd  that  the  flow  in  the  iiilei  or  onilei  port  reiii.nn- 
vijiioiiary  when  the  passage-end  eiicouiiiers  ihe  port  Foi  iln- 
reason  the  region  of  the  port  is  noi  included  in  the  s  -ni 
puiatioiial  dv<mam  shown  m  I  ig  2. 

The  time  dependent  process  ol  the  passage  end  iraiishnin* 
across  the  region  ivf  the  mici  or  otillei  port  will  be  relerri  )  .o 
the  "gradual  opening"  ol  the  pass.igo  The  pa-sage  opei.mg 
proeess  will  be  relerred  as  "iiisianianeous,"  when  ilu  .i- 
sumplioii  IS  made  that  the  pas-age  iiisiainaneou-ly  open-  o 
Ihe  port  area  and  is  siit'ieciej  iinmediaieh  lo  ihc  -leads  ilo.s 
conditions  at  ihe  port 

II  IS  assumed  ili.ii  the  flow  i-  mv  iseid  and  -an  be  modeh-il  'n 
Ihc  I  tiler  equations  Ihe  un-ieadv  iwo  dnnen-ion.il  I  ulei 
equations  can  be  w  riiien  m  coiiserv  ai  ion  law  lium  a- 
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Here  p  is  ihe  density,  u  and  v  the  velocity  components  in  the  X 
and  Y  coordinate  directions,  p  the  pressure,  and  >  the  ratio  of 
specific  heats.  7  he  energy  per  unit  of  volume,  e,  is  defined  by 
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where  t  =  P)iy- l)p  is  the  internal  energy.  We  look  for  the 
solution  of  the  system  of  equations  represented  by  Eq.  (1)  in 
the  computational  domain  shown  in  Fig.  2  in  time  /  with  the 
following  conditions  at  the  domain  boundaries:  a)  solid  wall 
along  segments  1-3  and  2-4,  b)  outflow  along  segment  3-4, 
and  c)  inlet  along  segment  1-2. 

It  is  assumed  that  initially  at  time  i  =  0,  the  passage  of  the 
wave  rotor  is  filled  with  stationary  gas  at  ambient  conditions. 
When  instantaneous  opening  of  the  wave  rotor  passage  was 
simulated,  it  was  assumed  that  at  time  f  =  0,  the  flow  at  the 
inlet  1-2  was  equal  to  the  steady  flow  in  the  port.  When 
gradual  opening  of  the  passage  was  simulated,  at  time  /  =  0, 
the  inlet  was  closed  and  solid  wall  boundary  conditions  were 
imposed  at  the  inlet  1-2.  Then,  this  boundary  condition  was 
gradually  replaced  by  the  flow  condition  at  the  inlet  port.  The 
length  uncovered  to  the  inlet  port  region,  where  the  solid  wall 
boundary  conditions  were  replaced  by  the  inlet  port  con¬ 
ditions,  was  determined  using  the  elapsed  time  and  the 
velocity  of  the  passage  relative  to  the  inlet  port. 

The  Godunov  method  was  used  to  obtain  a  numerical 
solution  of  Eq.  (1)  with  the  described  boundary  .ind  initial 
conditions.  Details  of  the  implementatioti  of  the  method  and 
boundary  conditions  are  given  in  Ref.  4. 

The  flowfield  was  modeled  for  the  rectangular  passage  with 
a  width  of  0.02  m  and  a  length  of  0. 12  m.  The  grid  cotering 
the  computational  domain  of  the  passage  is  shown  in  Fig.  2. 

Resull!)  and  Discussion 

The  following  initial  conditions  were  assumed  for  the  air 
initially  in  the  pas.sage: 

/’o=laim  po=1.2kg'm’  0,1  =  0  F<,  =  0 

The  driver  gas  entering  through  the  port  at  the  left-hand  end 
was  assumed  to  have  the  following  properties: 

P,,  =  l.8aim  Prf=l.8lkg/m’  L'rf=l50m^s  l',,=0 

These  condilions  correspond  to  a  practical  siluatiun  in  a  wave 
rotor  when  a  passage  filled  with  a  quiescent  fresh  charge  of  air 
at  ambient  condiiioiis  encounters  an  inlet  port  supplying  hoi, 
high-pressure  driver  gas. 

If  the  assumpiion  is  made  that  the  passage  insianiaiicously 
opens  and  is  subjected  to  the  condiiions  of  the  inlet  flow  ai  the 
Icl'i  boundary  1-2  (see  Fig.  2).  then  a  pericclly  onc- 
dimensional  flow  pattern  should  develop  in  the  passage.  I  hc 
results  of  the  modeling  of  these  couduions  are  presenied  m  the 
form  of  pressure  contours  at  progressively  larger  lime  steps  in 
Fig  3.  Time  /  =  0  corresponds  lo  the  moment  when  the 
passage  opens  liislaniancoiis  opening  of  the  passage  is  eeii  in 
Fig.  3  10  lead  lo  an  immediate  formation  of  a  shock  wave  and 
subsequent  propagaiinn  in  the  passage  from  ihe  lefi  lo  the 
right.  The  Mow  condiiions  at  the  inlet  port  are  matched  lo  the 
parameters  of  the  rarefaction  wave  which  clfcciively  will 


Fig.  2  Compiilaliiinal  diimain. 


T-0.023  msec  I 
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Fig.  3  Evolution  in  time  of  Ihr  shock  wivt  formed  after  inslaniane- 
ous  passage  opening. 


cancel  the  rarefaction  wave  moving  toward  the  inlet  at  the  left 
end  of  the  passage.  For  this  reason  the  flow  in  the  passage  has 
no  additional  disconiinuiiy.  This  situation  is  typical  for  a 
wave  rotor,  where  flow  condiiions  at  the  ports  are  chosen  in 
such  a  way  that  waves  do  not  propagate  from  the  passages 
into  inlet  or  outlet  ports.  The  flow  in  the  ports  will  iherefisre 
remain  steady.  In  the  case  shown  in  Fig.  3,  the  shock  wave 
was  found  lo  propagate  with  a  velocity,  f ,,,  =  in  s.  ilie 
interface  with  a  velocity  1  =  150  m  's,  and  all  the  paraineicrs 

esainnied  were  confirmed  accurately  using  one-dimensional 
gasdynamics  relationships.  , 

Most  of  the  approaches  used  in  wave  rotor  design  are  based 
on  the  assumption  lhai  the  waves  are  one-diineiiMoiial  in 
nature  Vkhen  the  gas  is  compressed  by  a  weak  sho^k  wave, 
lot  nisiance,  assninpiions  aie  made  ihai,  I)  the  pn'ccss  is 
isiKMiiropic,  2)  the  hot  and  cold  gases  are  strictly  sep.iraicd  by 
a  planai  intetface,  and  3l  the  flow  is  everywhere  irtv'taiional 
This  leads  to  the  very  high  efficiencies  projecicvl  lor  ihe 
compression  If  ihe  pass.igc  is  wide  enough  so  lhai  vissons 
cllecis  can  be  negicsied.  ihis  model  ol  ihe  compression  ni  ihc 
wave  rotor  passage  is  icahsiic.  bin  only  lor  insiaiii.incons 
passage  opening  or  lor  a  very  long  passage.  Resulis  prcsenicd 
below  illusiraie  how  ihe  gradual  passage  opening  aliccis  ihc 
compression  process 
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In  Fig.  4,  the  pressure  and  Maeli  number  contours  in  the 
passage  are  shown  at  a  sequence  of  times  for  the  case  ol  a 
gradual  opening  to  the  inlet  port  with  a  velocity  of  IIXJ  m/s. 
Parameters  for  the  gas  in  the  passage,  before  openiitg  begins, 
and  in  the  inlet  port  are  the  same  as  for  the  case  ol  in¬ 
stantaneous  opening.  The  dynamics  of  the  flow  dcvelopitient 
seeti  in  I  ig.  4a  is  very  different  from  that  shown  itt  Fig.  3. 
First,  curved  pressure  waves  radiaic  from  the  niitial  small 
opening  appearing  at  the  lower  corner  of  the  inlet  on  the  left 
side  of  the  passage  (/  =  0.044  ms  in  Fig.  4).  Subsequently, 
thes^waves  reflect  from  the  upper  wall  of  the  passage  and  at 
time  (  =  0.125  ms  have  formed  a  row  of  compression  waves 
which  have  approximately  straight  fronts  normal  to  the  wall 
of  the  passage.  Initially,  compression  waves  of  this  kind 
occupy  a  small  part  of  the  region  with  disturbed  gas.  In  the 
region  well  behind  the  front  of  this  quasi-one-dimensional 
propagation,  compression  waves  arc  curved  and  are  the  result 
of  the  interaction  between  the  waves  reflecting  back  and  forth 
between  the  lower  and  upper  walls  of  the  passage  and  the  new 
waves  created  by  the  progressive  opening  of  the  port  to  the 
passage.  Since  it  is  possible  to  see  in  Fig.  4  for  the  times 
(=0.084  and  0.166  ms,  respectively,  the  flow  behind  the 
quasi-one-dimensional  region  is  highly  rotational  and  is 
relatively  constant  pressure  in  the  ,V  direction.  The  passage 
became  fully  opened  only  at  (  =  0.2  ms.  At  this  time  com¬ 
pression  waves  are  propagating  along  the  length  of  the 
passage  and  the  pressure  rises  gradually  from  I  atm  at  the 
right  end  to  1.8  atm  at  the  left  end.  The  front  of  converging 
compression  waves  will  eventually  become  a  shock  wave,  but 
this  will  occur  at  some  later  time  and  outside  the  com¬ 
putational  domain. 

It  was  cottcludcd  frotn  the  case  presetited  iti  Fig.  4  that  a 
passage  leitgth-to-width  ratio  of  6  will  lead  to  very  high 
mixing  losses  and  noiiuniform  and  inefficient  compression, 
since  111  this  case  the  region  of  rotational  flow  occupies  half  of 
the  passage  volume.  Additional  losses  will  be  produced 
because  of  the  htghly  rotational  flow  withiti  each  of  the  two 
gases. 

In  Fig.  5.  pressure  and  velocity  contours  are  shown  for  ihe 
case  of  gradual  opening  of  the  passage  with  a  velocity  of  2tX) 
iti/s  l  ull  opening  of  the  passage  in  this  case  occurred  at 
(  =  0.1  Ills.  A  curved  shock  wave  is  formed  at  (  =  0.044  ms. 
This  shock  wave  (sec  Fig.  5a)  partially  reliccts  from  the  upper 
wall  ol  ilic  passage  and  ihcn,  gradually  converging  with  its 
mam  Irom.  forms  an  almost  planar  shock  (rout  at  the  lime 
(  =  0.252  ms.  However,  even  then  ihc  flow  is  highly  rotational 
behind  the  shock  from,  and  the  region  of  rotational  Mow 
occupis'v  one-third  of  the  passage  volume.  In  this  regisiii  the 
gas  velocity  increases  from  .5/  =  ().3  at  the  upper  wall  to 
M-  0.52  .11  the  lower  svall. 

W  lie  1  die  velociiy  ol  the  passage  opening  is  increased  to  5(X) 
in  s  (cec  I  le  6l  the  passage  becomes  fully  open  ai  /  =  0.04  ms. 
Because'  ol  the  last  opening  ot  the  passage,  the  shock  wave  at 
ilie  nine  '  =  0.04  ms  is  less  curved  and  only  a  small  tract  ion  ol 
ilie  shock  from  rellecis  Irom  (he  upper  wall.  Tins  reflected 
pan  ol  the  shock  Irom  converges  with  ilie  main  Irom  ai 
(  =  0  r2  ms.  From  that  lime  on,  the  How  paliern  in  the 
passage  is  mostly  one-dimensional  with  a  small  and 
weakening  region  ol  roiainnial  Mow  behind  the  main  Irom. 
.Ncveriheh'ss.  even  loi  this  case,  wnh  passage  lengrh  passage 
widili  .V  ihcrc  will  be  very  high  mivmg  losses,  wiili  a  large 
pan  ol  I  he  passage  volume  subjected  to  roiai tonal  How  . 

In  Older  lo  sindv  tu'w  the  sirengili  ol  Hie  shock  wave  m- 
llueiice'  'he  How  pailern  developing  in  Ihc  passage,  an  ad- 
diiinn.il  ,  ise  was  snmil.iied  where  Hie  paramclers  ol  Hie  dnver 
gas  ai  Ills  mici  port  area  were  mcreasecl;  namely,  to 

/’r=2.H5aini  t  ,  ~  2H3ms  p,  =  4kg'm' 

.\s  in  'he  pies  ions  e.ise.  Hie  paraineieis  of  ilic  driver  g.is  were 
chosen  in  siieli  a  w,iv  ih.ii  waves  dc>  noi  propagale  back  imo 
the  poll  shell  the  pass.ige  opens  Results  ol  line  sminhiiion 
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arc  shown  m  I  ig.  7.  The  velociiy  of  Hie  passage  opening  was 
I  =250  m  s.  It  was  concluded  from  the  results  presemed  in 
Fig.  7  that  an  increase  in  Hie  iiinial  shock  sireiigih  leads  lo 
stronger  rcHcctions  and  subsiatitial  increases  in  tlie  How 
roiaiion  I  Ins,  in  turn,  will  lead  lo  increased  losses  because  ol 
mixing  ol  the  driver  and  diiveii  gases.  The  pattern  ol  iiinlnplc 
relleciioiis  ol  the  shock  wave  Irom  the  nppci  and  lower  walls 
ol  the  passage  can  be  seen  in  I  ig  7a.  At  lime  (  =  0 O'f  nis  ilic 
shock  wave  rellccled  Irom  Hie  upper  wall  is  propagaimg 
ii'waul  ilie  lower  wall  I’arl  ol  Hic  rcflccied  shock  wave  liont 
converges  wnh  ilic  iiiaiii  shock  wave,  and  pan  begins  lo 
rclleel  Irom  the  lower  wall  ((  =  0.109  nis).  Ilie  iiuilliple 
refleciioii  weakens  the  rellecicd  shock,  but  ihe  retleciii’n 
belween  ihe  walls  ol  Hie  passage  cr'iilinues  and  can  be 
Ic'llowed  lor  (  '  0  2 1  1  and  0  245  ms.  Ii  is  clear  lhal  Un  iliesc 
ci'iuhliiuis  ecen  pass.iges  wnh  a  leiieili-lo.wiJih  raiio  ol  6  will 
have  snbstaniial  niixiiig  losses  because  ol  Hie  roiaiioii.d  I  low 
\  dll  leieiii  siinaiion  is  lonnd  lor  i  he  p.issage  opening  U'  I  he 
miei  pori  ol  Ihe  diiveii  gas  \l  ilns  pori,  ilie  pressure  and 
velociiy  ol  the  tdiicei)  gas  m  ihe  passage  should  be  maished 
wiih  Ihe  pressure  and  velociiv  of  ihc  (clioeii)  gas  ai  the  miei 
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walls.  Therefore,  the  interface  will  ''remeinber''  the  Jsnamics 
of  the  gradual  passage  opening  in  passages  wiih  large  lengih- 
lo-widih  ratio. 

Conclusions 

It  IS  shown  III  the  prescni  siujs,  on  the  bans  s>f  numeti..al 
modeling,  that  the  dsnamus  oi  me  passage  opening  sig- 
iiificanll>  affects  the  flow  paiieri,  i.  ile  passages  iit  wase 
rotor  desices.  Lor  rectangular  a>.j  passages,  esen  when  the 
velocity  of  the  passage  open  :  g  '•  -.  a  sine  dimtusional 
flow  pattern  forms  only  m  passage  wim  a  lengm  m-wijih 
ratio  larger  than  3  In  the  regi  n  pte.e,.!nig  ti-.  'otnia.ion  oi 
the  one-dimensional  fUiw  paiicm,  me  Ilow  .s  loiaiioiial  and, 
ill  some  instances,  contains  stunk  anj  pressures  waves  whuh 
repetitively  reflect  from  the  upper  and  lower  walls  ol  the 
passage.  In  practice  this  would  lead  lo  sigi  ilicani  mmng 
between  the  driver  gas  (e  g  ,  exhaust  gasi  and  driven  gas  le  g  , 
fresh  air)  and  reduce  the  elliciencj  ol  ihe  engine  svcte  W  iih  a 
reduction  in  the  velocity  of  the  passage  opening,  the  volume 
ol  the  mixing  region  increases  as  a  result  of  the  roiaiional 
flow  which  develops  from  the  slow  opening  of  ihe  passage 
With  an  increase  in  shock  wave  strength,  the  volume  of  ihe 
mixing  region  increases  as  a  result  of  the  rotational  flow 
which  develops  form  multiple  refleclions  of  ihc  shock  and 
pressure  waves  from  the  passage  walls. 

Numerical  modeling  of  the  process  of  gradual  passage 
opening  at  the  inlet  port  of  the  driven  gas  revealed  that  the 
interface  between  the  gases  will  move  all  the  way  through  the 
passage  with  a  frozen  pattern  of  distortion  or  obliqueness  lo 
the  passage  walls  The  interface  oblicyieiiess  increases  when 
the  velocity  of  (he  passage  opening  decreases. 


In  all,  it  can  be  concluded  that  taking  into  a^couiii  the 
gradual  passage  opening  is  essential  tor  the  wave  niactiinc 
design,  not  only  lor  proper  liming  and  wave  arrangeineni  but 
also  because  ol  the  losses  which  will  occur  due  to  inixing  and 
wave  rclleciioiis.  The  gap  between  projected  and  .ivhieved 
elficiencies  of  wave  machines-'''  could  be  pariiallv  due  lo 
neglect  of  the  effects  of  the  gradual  passage  opening  winch  are 
studied  herein. 
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The  problem  of  gradual  openina  of  rectangular  axial  passaae  in  wave 
rotors  was  studied  in  our  article^  published  in  the  January  issue  of  Journal 
of  Propulsion  and  Power.  There  we  have  defined  the  problem  of  qra<^”al 
opening  in  the  wave  rotor  passaae  and  the  mathematical  model  which  was  used 
to  simulate  the  opening  process.  In  this  article  we  will  assume  that  the 
reader  is  familiar  with  the  reference^,  and  it  could  be  reoarded  as  an 
extension  to  the  that  reference. 

If  the  wave  rotor  is  used  to  produce  shaft  power,  it's  passaoes  should 
be  skewed  in  one  form  or  another.  In  this  paper  we  will  analyze  the  aradual 
openina  of  a  skewed  passaae  and  will  examine  the  conclusion  drawn  for  the 
rectanoular  axial  passaae  for  this  more  aeneral  aeometrv  of  the  passaae. 

{Results’  and'  biscussibh 

The  main  conclusion  of  the  study  on  oradual  openina  of  rectanoular 
passaae  is  that  in  order  to  minimize  the  mixing  losses  caused  by  rotational 
flow  in  the  passaoe,  the  openino  velocity  should  be  very  hioh.  In  the 


limit,  instantaneous  openina  of  the  wave  rotor  passaoe  will  lead  to  one 


dimentional  flow  pattern  in  the  passaae  will  have  minimal  mixina  losses. 

When  the  passaae  of  the  wave  rotor  is  skewed,  even  instantaneous  openina  of 
the  passaae  will'  hb^  lead  to  development  of  the  one  dimensional  flow  pattern 
with  small  mixina  losses.  That  statement  will  be  demonstrated  in  the 
followina  example. 

Let's  model  the  openino  process  for  the  passaoe  0.02  m  wide  and  0.24  m 
Iona.  It  has  left  and  riaht  hand  inlets  parallel  to  the  v-axis  and  upper 
and  lower  walls  of  the  passaae  form  60®  anale  with  positive  direction  of  the 
x-axis.  It  is  assumed  that  initially  air  in  the  passaae  is  at  the  followina 
conditions . 

P  ■  1  atm;  P  *1.2  ka/M^ ;  -  0;  V  *  0 

O  o  '  0  O 

The  driver  oas  enterina  throuoh  the  port  at  the  left  hand  end  was  assumed  to 
have  the  followino  properties; 

P.  *  1.P  atm;  P.  *  1.P1  ko/M^;  »  75  M/sec;  =  129.9  M/sec 
d  .  d  d  d 

The  conditions  for  the  driver  and  driven  aas  are  the  same  as  in  case  of  the 
rectanoular  passaae  which  we  reported  in  the  previous  paraoraph. 

In  Fioures  la  and  1b  results  for  simulation  of  the  instantaneous 
openina  of  the  skewed  passaoe  are  sh<>m  in  form  of  pressure  and  velocity 
contours  at  the  seauences  of  times.  The  flow  pattern  near  the  inlet  in 
Fiaure  1  is  hiohly  rotational  which  suooests  very  hioh  mixina  losses.  That 
is  caused  partially  by  reorientation  of  the  shock  wave.  At  the  first 


moments  of  time  after  the  openinq  of  the  passaqe  the  shock  wave  between  the 
driven  and  driver  oases  is  oblique  to  the  lower  and  upper  walls  of  the 
passaqe.  At  later  times  this  shock  wave  turns  and  becomes  normal  to  the 
upper  and  lower  walls.  Thus,  incontrast  with  rectanqular  for  skewed 
qeometry  passaoes  there  is  no  obvious  condition  for  minimizino  the  mixino 
losses  caused  by  the  inlet  openinq. 

Let's  find  the  conditions  for  openinq  of  the  skewed  passaqe  which  will 
lead  to  minimal  mixino  between  driver  and  driven  qases.  As  we  have  stated 
before,  rotational  flow  in  the  instantaneously  opened  skewed  passaqe,  will 
develop  because  of  the  rotation  of  the  shock  wave  from  the  parallel  position 
to  the  inlet  (initially)  to  the  position  where  it  is  normal  to  the  upper  and 
lower  walls  of  the  passaqe  (when  the  flow  is  developed).  If  we  could  form  a 
shock  wave  which  will  be  all  the  time  normal  to  the  lower  and  upper  passaoe 
walls  then  the  rotational  flow  and  mixinq  will  be  minimal. 

Analysinq  the  formation  of  the  shock  wave  front  at  the  inlet  we 
concluded  that  for  the  shock  wave  formino  at  the  inlet  to  remain  normal  to 
the  lower  wall,  inlet  should  be  openinq  at  the  rate: 

%  ■  ''sh  > 

where  V  -  velocity  of  the  inlet  openinq. 
op 

-  shock  wave  velocity  in  the  media. 

-  anqle  of  the  skewed  passaqe . 

The  velocity  is  equal  to  the  velocity  with  which  the  shock  wave  surface  will 
slide  alonq  the  inlets  wall. 

In  Fiqs .  2a  and  2b  results  for  the  modelinq  of  the  oradual  openinq  of 
the  skewed  passaqe  inlet  with  the  openinq  velocity  in  accordance  to  equation 


(2)  are  shown.  As  it  is  obvious  from  coutour  plots  of  pressure  and  velocity 

in  Fiqs.  2a  and  2b,  openinq  the  passaoe  with  velocity  V  calculated  from 

op 

equation  (1)  assures  development  of  the  flow  pattern  with  minimal  mixina. 

The  shock  wave  surface  in  that  case  remains  at  all  times  normal  to  the  walls 
of  the  passaqe. 

bohcliisiohs 

Modelinq  of  the  qradual  openinq  of  the  skewed  passaqe  revealed  the  way 
to  minimize  mixinq  losses  at  the  passaae  inlet.  The  mixinq  will  be  minimal 
when  the  velocity  of  the  openinq  is  matched  with  shock  rave  velocity  in  the 
passaqe  divided  by  the  anqle  of  skewed  passaqe.  Our  simulations  shews  that 
even  v^en  conditions  for  the  optimal  openinq  are  not  satisfied  exactly, 
reduction  of  the  rotational  mixinq  could  be  siqnificant  if  the  openinq 
velocity  is  jf15%  of  the  optimal  value. 

We  can  also  conclude  that  a  very  hiqh  openinq  velocity  is  required  in 
order  to  obtain  low  mixinq  loss.  The  openinq  velocity  of  the  passaqe  should 
be  always  hioher  than  the  velocity  of  the  shock  wave  qenerated  in  the  wave 
rotor  passaqe  at  the  hiqh  pressure  inlet.  Fven  substantially  skewed  passaae 
will  require  the  openinq  speed  *10%  hioher  than  the  shock  wave  velocity  in 
the  passaqe  which  is  not  easy  to  achieye  for  the  typical  flow  conditions  in 
the  wave  rotor. 

Another  limitation  is  that  even  when  this  optimal  speed  of  openino  is 
achieved  at  one  port,  it  will  not  be  optimal  at  other  ports;  so, 
optimization  will  not  be  full.  However,  mixinq  losses  in  the  skewed 
passaqes  will  be  smaller  than  in  rectanoular  one,  if  the  qradual  openino  of 
the  skewed  passaae  beqins  from  its  acute  anole. 
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ABSTRACT 

The  dynamics  of  the  expansion  of  detonation 
products  la  a  diverging  channel  were  Investigated 
numerically.  A  two-dimensional  unsteady  Euler  code 
based  on  the  Godunov  method  was  employed.  The 
Influences  of  particular  features  In  the  expansion 
process,  such  as  the  presence  of  reversed  flow  and 
propagation  of  hammer  shocks,  on  the  production  of 
thrust  were  examined.  Sequentially  expansions  of 
detonation  products  were  also  modeled  and  it  was 
concluded  that  In  order  to  maintain  a  high  frequency 
periodic  mode  of  operation  for  propulsion  applications 
the  channel  should  be  refilled  with  ambient  air  after 
each  expansion.  The  Influence  of  the  ratio  of  ambient 
air  to  detonation  product  volume  on  the  dynamics  of  the 
thrust  production,  and  on  the  Impulse  generated  during 
the  expansion,  are  also  reported. 

NOMENCLATURE 

dAjj  -  element  of  channel  Internal  surface  area  normal 
to  the  x-axls 

e  -  energy  per  unit  volume 

p  “  pressure 

t  -  time 

u  -  component  of  velocity  In  x  direction 

V  -  component  of  velocity  in  y  direction 

X  -  cartesian  coordinate 

y  “  cartesian  coordinate 

e  -  Internal  energy  per  unit  mass 

y  -  ratio  of  specific  heats 

p  -  density 

subscripts 

0  -  conditions  at  t-0  in  combustion  products 

ex  -  external 

L  -  along  Inner  wall  of  the  channel 
INTRODUCTION 

Early  In  the  development  of  propulsion  engines  for 
aircraft,  a  choice  had  to  be  made  between  engine 
concepts  based  on  nonsteady  or  on  steady  gasdynamlc 
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processes.  Concepts  based  on  steady  gasdynamlc 
processes  were  pursued  largely  because  they  were 
simpler  to  analyze  and  to  appreciate  conceptually.  On 
the  basis  of  the  chosen  steady  concepts,  engines  were 
built,  developed  and  perfected  with  time.  The 
nonsteady  concepts,  which  looked  as  promising 
Initially,  were  not  developed  and  remain  today  In  the 
conceptual  stage. 

The  need  for  a  more  thermally  efficient  small 
engines  for  various  military  and  civilian 
applications,  has  led,  at  Intervals,  to  a 
reexamination  of  nonsteady  concepts  (1 ,2,3).  Today, 
such  an  examination  can  be  made  more  thoroughly.  For, 
while  nonsteady  engine  concepts  were  at  one  time 
extremely  difficult  to  analyze,  numerical  modeling  on 
computers  can  now  provide  time-dependent  pictures  of 
Internal  flow  processes.  This  can  provide  efficient 
tools  for  preliminary  design  calculations  and 
eventually  analysis  tools  for  design  optimization. 

From  the  late  50's  until  the  early  70's,  the 
feasibility  of  an  engine  operating  with  Intermittent 
detonatlve  combustion  was  studied  at  the  University  of 
iMlchigan  (A, 5, 6).  Specific  Impulses  over  2100  sec, 
were  realized  for  a  single  linear  shock  tube  operating 
Intermittently  with  frequencies  up  to  35  detonations 
per  second.  Nicholls  et  al.  showed  that  an  engine 
operating  Intermittently  on  detonation  waves  will  have 
some  advantages  over  an  engine  operating  on 
deflagration.  In  general,  it  will  have  very  high 
specific  thrust.  The  engine  will  be  very  simple 
mechanically,  and  not  require  precompression  of  the 
mixture  for  efficient  combustion. 

One  of  the  disadvantages  of  the  engine  Is  that  Jet 
velocities  developed  after  detonatlve  combustion  are 
very  high.  Efficient  operation  of  the  engine  for 
propulsion  could  be  realized  only  at  high  supersonic 
vehicle  velocities  (Mach  number  of  about  4).  In  order 
to  reduce  the  velocity  of  the  out-coralng  Jet,  it  was 
proposed  to  use  spinning  detonation  (^).  But  the 
spinning  detonation  process  proved  to  be  unsuitable 
for  use  In  an  engine  because  It  was  unstable.  Back  et 
al.  (^)  studied  the  feasibility  of  using  detonatlve 


propulsion  In  Jupiter  s  atmosphere,  and  obtained  a 
number  of  results  for  various  types  of  nozzles. 

Motivated  by  the  potential  applications,  of 
detonatlve  propulsion,  In  the  present  study  the 
process  of  the  expansion  of  detonation  products  In  a 
diverging  channel  or  diffuser  of  a  detonation  engine 
were  modeled  numerically.  The  model  and  computer 
program  were  then  applied  to  calculate  repetitive 
firing,  and  the  effect  of  changing  the  diffuser 
length. 

THE  MATHEMATICAL  MODEL  AND  NUMERICAL  SOLUTION 

It  Is  assumed  that  the  jet-propulsion  nozzles  of 
an  engine  using  detonatlve  combustion  can  be 
constructed  from  a  number  of  straight  diverging 
channels  which  are  grouped,  or  clustered  together.  In 
one  cycle  of  the  engine,  a  small  volume  of  each  tube 
la  filled  with  the  combustible  mixture  and  undergoes 
detonation.  Then  the  detonation  products  expand  into 
the  stationary  gas  which  fills  the  rest  of  the  tube  at 
ambient  pressure  and  temperature.  The  tube  is 
diverging  and  so  the  shock  wave  decays  in  strength. 
Finally,  a  weak  shock  wave  leaves  the  tube,  and  a 
subsonic  flow  of  gas  discharges  from  the  exit  of  the 
channel. 

It  is  assumed  that  the  flow  Is  periodic  along  the 
lines  which  are  a  continuation  of  the  walls  of  the 
channel.  This  implies  that  we  are  modeling  the  flow 
In  a  centrally  located  channel  of  the  cluster. 

It  Is  further  assumed  that  the  flow  Is  Invlscld. 

The  unsteady  two-dimensional  Euler  equations  can 
be  written  In  conservation  law  form  as: 


1£ +  1I  +  1£-  0 

at  ax  aY 


where 
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p  +  pv^ 
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Here  p  Is  the  density,  u  and  v  are  the  velocity 
components  In  the  X  and  Y  coordinate  directions,  p  is 
the  pressure  and  y  Is  the  ratio  of  specific  heats. 

The  energy  per  unit  of  volume,  e.  Is  defined  as 

2.2 

e  -  p(  e  + - 2 - )  • 

where  t  ■  ^  is  the  Internal  energy  per  unit  mass, 

We  look  for  the  solution  of  the  system  of  equations 
represented  by  Eq  (1)  in  the  computational  domain  as 
shown  In  Fig  1  in  time  t,  with  the  following 
conditions  at  the  domain  boundaries: 


a)  Solid  wall  along  segment  1-2,  1-3  and  2-4 

The  condition  on  the  surface  is  defined  by  solving 
the  Rlemann  problem  between  the  point  nearest  to 
the  wall  in  the  domain  of  integration  and  Its 
mirror  image  in  the  direction  normal  to  the  wall, 

b)  Periodicity  between  segments  3-5  and  4-6 

By  virtue  of  the  central  location  of  the  channel 
flow  at  each  grid  point  along  3-5  is  the  same  as 
at  corresponding  grid  points  along  4-6  in  Fig  1. 


c)  Outflow  along  segment  5-6 

If  the  outflow  is  subsonic  at  the  downstream 
boundary  we  define  only  the  pressure  pQut* 
for  all  other  flow  parameters  apply  the 
continuation  condition.  That  means  that  Uq^^, 
Vout  Pout  equal  to  the  values  of  u,  v 

and  p  one  point  ahead  of  the  downstream 
boundary. 

If  the  outflow  is  supersonic  the  continuation 
condition  is  applied  to  all  parameters  at  the 
downstream  boundary. 

It  is  assumed  that  initially  at  the  time  t-0,  the 
channel  is  filled  with  detonation  products  from 
segment  1-2  to  l'-2'.  Outside  this  region  of  the 
computational  domain,  the  gas  is  assumed  to  be  at 
ambient  conditions. 

The  Godunov  method  has  been  used  to  obtain  a 
numerical  solution  of  the  Eq  (1)  with  the 
described  boundary  and  Initial  conditions.  Details  of 
the  Implementation  of  the  method  and  boundary 
conditions  are  given  in  (^). 

An  orthogonal  grid  is  constructed  which  covers 
the  computational  domain  as  shown  on  Fig  1.  It  should 
be  noted  that  although  the  formulation  of  the  problem 
and  its  solution  is  two-dimensional,  with  the  boundary 
conditions  described  above  the  flow  will  be 
essentially  one  dimensional. 

It  is  assumed  that  at  t"0,  the  detonation  wave 
has  passed  through  the  combustible  mixture  and  the 
products  of  combustion  have  uniform  properties.  Thus 
the  time  of  initiation  and  propagation  of  the 
detonation  wave  through  the  combustible  mixture  is 
neglected  In  comparison  with  the  time  for  propagation 
through  the  channel  and  subsequent  expansion  of  the 
gases. 


RESULTS  AND  DISCUSSION 


The  history  of  the  flow  development  in  the 
diverging  channel  is  shown  in  Fig  2.  Here  the 
velocity  and  pressure  distributions  at  the  center-line 
of  the  channel  are  shown  at  specific  moments  of  time, 
where  time  t»0  is  when  the  expansion  begins.  At  time 
t»0,  approximately  8.5X  of  the  channel  volume  from  the 
left-hand  side,  is  filled  with  detonation  products 
with  initial  pressure  density  and  velocity  having 
values  of 

Po  “  40  atm,  Pq  ■  5.2kg/m^,  Vq  ■  0  and  Uq  ■  0 
respectively.  At  t*0,  the  rest  of  the  channel  is 
filled  with  ambient  air  with 

p  ■  1  atm,  p  *  1.3  ^8/m^,  v  ■  0  and  u  ■  0  , 

The  length  of  the  channel  la  I  meter.  These 
conditions  will  be  referred  to  as  Case  1. 

Initially,  a  shock  wave  going  from  left  to  right 
and  a  rarefraction  wave  going  from  right  to  left  can 
be  seen  in  Fig  2.  The  shock  wave  progressively  decays 
because  the  channel  is  diverging.  The  rarefraction 
wave  propagates  rapidly  towards  the  wall  at  the  left 
end  of  the  channel,  because  the  temperatures  in  this 
region  are  very  high  (approximately  2800  K)  and  the 
channel  is  converging  towards  the  left  end.  At  a  time 
of  about  0,15  msec  the  rarefraction  reaches  the  wall 
at  the  left  end  of  the  channel.  From  this  moment, 
pressure  at  the  left  end  of  the  channel  rapidly 
decreases  from  about  40  atm  at  t  ■  0.15  msec  to  about 
0.3  atm  at  t  ■  0.82  msec.  Because  of  this  rapid 
depressurization  of  the  region  adjacent  to  the  left 
end  wall,  the  pressure  on  the  right  hand  end  of  the 
channel  becomes  higher  than  that  at  the  left,  which 
generates  a  shock  wave  going  from  right  to  left  at 

time  toO.6  msec.  This  back-going  shock  wave  rapidly 
slows  down  and  then  reverses  the  direction  of  the  flow 
because  the  channel  is  converging  towards  its  left  end 
wall. 

The  negative  flow  reflects  from  the  left  end  wall 
at  time  t*<1.4  msec,  forming  a  hammer  shock  wave  which 
moves  from  left  to  right  and  again  reverses  the  flow 
direction.  The  formation  of  the  hammer  shock  wave  can 
be  seen  clearly  on  the  velocity  and  pressure  graphs 
for  time  t  ■  1,6  msec.  Because  the  shock  travels  in  a 
diverging  channel,  it  weakens  rapidly. 

The  main  shock  wave  produced  by  the  detonation 
products  leaves  the  channel  at  the  time  ■0.94  msec, 
and  because  the  channel  is  diverging  in  the  direction 
of  propagation  the  pressure  behind  the  shock  wave 
drops  from  about  10  atm  at  time  t"0.1  msec  to  about  5 
atm,  at  t  ■  0.94  msec  when  it  leaves  the  channel. 


The  mass  velocity  of  the  gas  drops  from  about 
1100  m/sec  at  the  beginning  of  the  expansion  to  about 
500  m/sec  when  the  gas  begins  to  leave  the  detonation 
tube.  The  velocity  of  the  gas  leaving  the  tube  then 
decreases  Co  -20  m/sec  and  Increases  again  to  about 
50  m/sec  when  the  shock  wave  reflected  from  the  left 
end  wall  reaches  Che  exit  of  Che  channel.  This 
concludes  the  significant  events  which  occur  in  the 
diverging  channel  as  a  result  of  Che  expansion  of  the 
detonation  products.  At  later  times,  rarefractlon  and 
pressure  waves  continue  to  form  and  travel  in  the 
channel,  but  they  are  much  weaker  and  their  influence 
on  the  thrust  produced  is  very  small. 

In  Fig  3  thrust  as  function  of  time  is  plotted 
for  the  case  described  in  the  graphs  given  in  Fig  2. 
Thrust  was  calculated  for  the  case  of  a  stationary 
tube  with  external  pressure  equal  to  the  ambient 
pressure  using 

"  -  -  Pex>^x  (2) 

It  can  be  seen  in  Fig  3  that  most  of  the  thrust 
is  produced  in  the  2  msec  following  the  beginning  of 
the  expansion  of  the  detonation  products.  Integration 
of  thrust  in  time  shows  that  902  of  the  impulse  is 
produced  in  the  first  2  msec  and  102  in  the  following 
7  msec. 

The  peculiarities  of  the  thrust  curve  can  be 
understood  in  relation  to  the  wave  pattern  in  the 
tube.  The  first  change  in  slope  at  time  t-O.S  msec 
corresponds  to  the  tine  when  the  front  of  the  main 
shock  wave  reaches  the  end  of  the  tube.  Because  the 
pressure  behind  the  shock  front  drops  rapidly,  the 
thrust  decay  increases  when  the  main  shock  front 
leaves  the  channel.  At  time  t“1.5  msec  the  thrust 
increases  as  a  result  of  the  hammer  shock  wave 
reflecting  from  the  left  wall. 

When  the  main  process  of  expansion  ends,  at 
time  t"6  msec,  the  pressure  in  the  tube  is  about  l.l 
atm  and  the  average  density  is  about  0.43  kg/m^.  This 
corresponds  to  an  average  temperature  in  the  channel 
of  about  935  K.  Since  the  pressure  in  the  tube  is 
close  to  ambient,  the  heat  exchange  will  be  by  natural 
convection,  and  it  will  take  a  relatively  long  time 
for  the  gas  remaining  in  the  tube  after  the  expansion 
to  cool. 

In  order  to  examine  what  will  be  the  wave  pattern 
and  thrust  when  a  second  detonation  wave  expands  in 
the  channel  immediately  following  the  expansion  of  the 
first  detonation  wave  (6  msec  after  the  beginning  of 
the  process  of  expansion  of  the  first  detonation 
wave),  conditions  were  simulated  numerically  in  the 
following  way.  Results  of  the  calculation  for  the 
first  detonation  wave  expansion  were  used  as  the 
initial  distributions  of  flow  parameters.  Then,  in 


the  8.5Z  of  the  channel  volume  adjacent  the  left  end 
wall  the  pressure,  density  and  velocity  were  given  the 
values 

Pq  ■  40  atm,  -  5.2  and  Vq  ■  Uq  ■  0 

respectively  at  time  t*0.  These  conditions  will  be 
referred  to  as  Case  2. 

In  Fig  4,  the  thrust  vs  time  Is  shown  for  the 
second  detonation  wave  expansion.  Since  the 
detonation  products  expand  In  a  low  density  hot  gas, 
the  expansion  goes  much  faster  than  In  the  previous 
case.  At  time  t*0.8  msec  from  the  beginning  of  the 
second  detonation  wave  expansion,  the  thrust  drops  to 
zero.  This  can  be  compared  with  a  thrust  of 
approximately  6*10^N/m  at  the  corresponding  time  shown 
in  Fig  3,  for  expansion  Into  ambient  air.  Then, 
because  the  large  volume  of  the  channel  has  static 
pressure  lower  than  ambient ,  thrust  becomes  negative 
and,  as  In  the  previous  case,  a  back-going  shock  wave 
develops.  In  Fig  5,  the  development  in  time  of  this 
recursive  shock  can  be  followed.  It  Is  observed  that 
the  negative  flow  velocities  which  develop  are  twice 
as  large  as  Chose  for  the  preceding  expansion  Into 
ambient  air.  And  at  time  t>>2  msec  the  flow  throughout 
Che  tube  Is  reversed.  Reflection  of  Che 
left-traveling  shock  wave  from  the  left  end  wall  at 
time  tx2  msec  Che  flow  throughout  the  Cube  la 
reversed.  Reflection  of  the  left-traveling  shock 
wave  from  the  left  end  wall  at  time  to2  msec  produces 
a  very  strong  hammer  shock  which  reverses  the  thrust 
and  the  flow  direction  at  time  to2.1  msec.  Because  of 
the  large  negative  thrust,  the  total  Impulse  which  Is 
produced  In  Case  2  Is  about  15Z  smaller  than  the 
Impulse  of  the  detonation  products  expanding  In 
ambient  air. 

Another  notable  effect  of  expansion  Into  the 
products  of  Che  previous  cycle.  Is  the  substantial 
rise  of  temperature  which  occurs  In  the  channel  after 
the  second  expansion.  When  Che  process  of  the  second 
expansion  ends  the  average  temperature  In  the  channel 
Is  approximately  1550  K.  This  Is  615  K  higher  than 
Che  final  temperature  after  expansion  into  ambient 
air. 

In  order  to  examine  how  thrust  and  total  Impulse 
are  Influenced  by  Che  volume  of  air  contained  at 
ambient  conditions  In  Che  channel  before  Che 
expansion,  two  additional  test  cases  were  modeled.  In 
Case  3,  the  volume  of  ambient  air  in  the  channel  to 
the  right  side  of  the  interface  with  the  detonation 
products  is  taken  Co  be  approximately  2.5  times  larger 
than  the  volume  of  the  detonation  products.  This 
requires  a  channel  about  50Z  shorter  than  the  original 
channel,  for  which  the  ratio  of  the  volumes  of  ambient 
air  to  detonation  products  was  about  11.  In  Case  4, 
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there  was  no  air  in  Che  channel  when  expansion  began, 
which  Implies  that  the  channel  In  this  case  has 
only  about  25Z  of  the  original  length.  Therefore, 
between  Cases  3  and  4  and  Case  1  only  Che  length  of 
the  channel  was  changed.  The  geometry  of  the  channel. 
Initial  parameters  of  the  detonation  products  were  the 
same  In  Cases  1,  3  and  4. 

In  Fig  6  the  thrust  Is  shown  as  a  function  of 
time  for  Cases  1,  3  and  4.  It  can  be  seen  that  Case  1 
(solid  line)  where  Che  channel  Is  fully  extended, 
gives  consistently  higher  thrust  and  that  the  thrust 
decreases  when  the  volume  of  ambient  air  In  the 
channel  (or  channel  length)  decreases.  Case  3  (chain- 
dotted  line)  and  Case  4  (dashed  line)  have  regions 
where  thrust  Is  negative.  Comparison  of  the  total 
Impulse  generated  shows  Chat  the  extension  of  Che 
channel  Is  very  beneficial  for  producing  Impulse.  The 

Impulse  produced  In  Case  4,  where  Che  Cube  Is  filled 
with  detonation  products  la  0.162*l(pN/M*8ec.  When 
the  channel  Is  extended  so  that  only  28.6%  of  the 
channel  volume  Is  Initially  filled  with  detonation 
products  (Case  3),  Che  Impulse  Increases  28%.  An 
additional  extension  of  the  channel  so  that  only  about 
8.5%  of  the  tube  volume  Is  Initially  filled  with 
detonation  products  (Case  1),  leads  to  total  impulse 
of  0,2784»103  N/M’sec,  which  is  a  72%  increase 
over  Case  4. 

CONCLUSIONS 

Numerical  modeling  of  the  process  of  expansion  of 
detonation  products  In  a  diverging  channel  revealed  a 
flow  pattern  rapidly  changing  in  time  with  multiple 
shock  and  rarefractlon  waves.  One  of  the  most 
Interesting  features  was  Che  occurrence  of  a 
backwards-traveling  shock  wave,  which  developed  as  a 
result  of  over-expansion  In  the  channel.  The  shock 
wave  first  reversed  the  flow  direction  (and  In  some 
cases  developed  a  negative  thrust)  and  then, 
reflecting  from  the  left  end  wall  as  a  hammer  shock, 
reversed  the  flow  direction  again.  Increasing  the 
thrust. 

Modeling  of  the  second  detonation  wave  expanding 
Into  the  products  of  the  previous  expansion,  showed 
that  this  arrangement  leads  to:  a)  significant 
negative  thrust;  b)  Increases  of  the  gas  temperature 
after  expansion;  c)  losses  of  15%  of  the  total 
Impulse.  This  leads  to  the  conclusion  that  the 
channel  should  be  filled  with  ambient  air  after  each 
detonation  and  expansion.  Since  the  process  of 
expansion  takes  only  about  0.006  sec,  the 
Clme-averaged  thrust  which  can  be  generated  by  a 


single  channel  will  be  limited  by  the  rate  at  which 
the  channel  can  be  re-charged  with  air  and  mixture, 
and  reflred.  For  very  low  Initial  velocities  In  the 
channel  (consistent  with  the  Initial  conditions 
specified  In  the  modelling),  the  frequency  of  firing 
Is  limited  to  about  20  detonations  per  second.  This 
corresponds  to  0.044  sec  for  recharging  at  gas 
velocities  of  about  23  m/sec.  Investigation  of  the 
Influence  of  the  channel  length  on  the  Impulse  and 
thrust  production  showed  that  Increases  In  the  volume 
of  ambient  air  In  the  channel  ahead  of  the  detonation 
products,  which  was  obtained  with  an  Increase  In 
channel  length,  led  to  significant  Increases  In  the 
thrust  and  Impulse  generated  during  the  expansion. 

The  total  Impulse  generated  by  an  expansion  In  the 
channel  In  which  8.52  of  the  total  volume  was  filled 
Initially  by  the  detonation  products,  was  found  to  be 
722  greater  than  In  the  case  of  a  short  channel  filled 
only  with  the  same  volume  of  detonation  products  at 
the  same  initial  conditions. 

The  diverging  channel  geometry  provided  the 
simplest  geometry  with  area  change  to  model  using  the 
unsteady  2D  Godunov-Euler  code.  Variations  In  the 
geometry  can  now  be  made  In  order  to  examine  the 
unsteady  thrust  performance. 
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FIGURES 


Figure  4.  The  Thrust  vs.  Time  Curve  for  the  Second 
Sequential  Expansion  of  the  Detonation 
Products. 
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INTRODUCTION 

Wave  Rotors  are  devices  which  use  unsteady 
pressure  waves  to  effect  direct  energy  exchange  between 
two  gates.  They  offer  the  potential  for  significant 
laproveoents  In  alr-'breathlng  engine  propulsion  cycles 
through  their  capability  of  withstanding  higher 
teaperatures  and  pressures  than' present-day  turbines. 
Other  diverse  applications  are  recounted  In  Ref  (1). 

In  Its  aost  basic  fora,  a  wave  rotor  It  a  drua 
with  axial  or  haltcal  (staggered)  passages  arranged 
around  the  periphery.  This  tingle  drua-llke  rotor 
replaces  separate  coapressor  and  turbine  coaponentt  for 
gas  turbine  applications.  The  coaprcsslon  and 
expanalon  of  the  two  fluids  occurs  In  the  passages  ss  a 
result  of  shock  tube  like  processes.  In  a  typical 
configuration,  coabustlon  products  at  high  teaperature 
and  pressure  give  up  energy  to  the  'driven'  fluid 
(usually  air)  through  the  action  of  tlae-unsteady 
coapresslon/thock  waves.  The  coabustlon  products,  (hot 
'driver'  gates),  are  In  turn  expanded  and  exhausted 
froa  the  rotor,  the  available  work  of  expansion  being 
utilised  to  Induce  a  fresh  charge  of  air  onto  the 
rotor.  Careful  sequencing  of  the  passage  ends  past 
stationary  Inlet  and  outlet  ports  In  valve  plates  on 
either  side  of  the  rotor  creates  a  cyclic  Internal  wave 
pattern  In  the  wave  rotor  coaponent.  The  high 
teaperature  capeblllty  of  Che  device  Is  a  direct 
consequence  of  the  repetitive  processing  of  both  cold 
and  hot  fluid  alternately  In  the  saae  rotor. 

Estlaatlon  of  the  aero-theraodynaalc  perforaance  of 
these  devices  hinges  on  calculations  of  Che  unsteady 
energy  exchange  In  the  wave  rotor  coaponent.  Nuaerlcal 
slaulatlon  of  the  unsteady  wave  processes  can  generally 
be  carried  out  on  a  one-dlaenalonal  basis.  However, 
the  coaplex  pattern  of  flow  discontinuities  and  wave 
Interactions  known  to  exist  In  wave  rotors  call  for 
nuaerlcal  asthods  which  solve  the  nonlinear,  hyperbolic 
systea  of  governing  equations  without  relying  on  cither 
artificial  viscosity  or  special  treacaenc  of 
discontinuities.  Described  In  this  paper  are  two  such 
techniques,  the  Godunov  aechod  and  a  Randoa  Choice 
aethod  (Cllaa's  aechod  or  Piecewise  Saapllng  aechod). 


NUMERICAL  PORMULATION 

An  essential  coaponent  of  both  the  Randoa  Choice 
aethod  and  the  Godunov  aethod  Is  the  solution  of  a 
aequence  of  ’;>cal  Rleaann  probleas.  A  Rleaann  problea 
la  defined  by  setting  up  an  Initial  value  problem  for 
the  equations  of  aotlon  In  Eulerlan  fora.  The 
unsteady,  cwo-dlaenslonal  Euler  equations  can  be 
written  In  vectorised,  conservation  law  fora  as: 

Ut  +  lP(U))x  +  lC(U))y  •  0,  where  U  •  1  p  (u  pv  e)'^, 

F(U)  “  [  pu  (p  +  pu^)  puv  (e  *  p)u)^,  and 

C(U)  ••  Ipv  puv  (p  +  pv^)  (e  +  p)vlT  . 

Here,  p  Is  Che  density,  u  and  v  are  Che  velocity  coa- 
ponents  In  Che  x  and  y  coordinate  directions,  and  p  Is 
Che  pressure.  Subscripts  Indicate  partial  differenti¬ 
ation  with  respect  to  that  Independent  variable.  The 
energy  per  unit  volume,  e,  Is  defined  by: 

e  ■  p(  c  (u^  v^)/2),  where  c  •  p/{'T  -l)p 

Is  the  Internal  energy  per  unit  aass. 

The  Rleaann  problea  Is  Intrinsically 
one-dlaenslonal  In  nature  and  Is  defined  accordingly. 
For  Che  one  space  variable  case,  thus,  U  -  U(X,c),  and 
Che  Initial  value  problem  Is  set  up  by  specifying 
Initial  conditions  which  consist  of  Intervals  over 
which  data  are  constant,  separated  by  jump 
discontinuities.  If  Che  tlae  step  Is  chosen  to  be 
sufficiently  small,  Che  waves  propagating  with  finite 
speeds  froa  adjacent  discontinuities  remain  within 
their  respective  spatial  cells  and  do  not  Intersect. 
This  sequence  of  solutions  froa  adjacent  Rleaann 
probleas  Is  then  pieced  together  to  obtain  Che  whole 
solution  at  each  succeeding  tlae  step. 

The  Randoa  Choice  aechod  (RCM)  for  gaadynaalcs  la 
the  outcoae  of  a  constructive  existence  proof  of 
solutions  Co  systems  of  nonlinear  hyperbolic 
conservation  laws  presented  by  Gllns  Ref  (2),  and  Its 
developaent  Into  an  effective  nuaerlcal  tool  by  Chorln 
Refs  (3,4).  The  solution  Is  advanced  In  time  by  a 
aethod  Chat  Includes  a  solution  of  Rleaann  probleas  as 
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datcrlbcd  above  and  a  aaapllng  procedure,  which 
chooaea  valuea  frow  a  repreaentaclve  point  of  the 
exact  aolutlon  to  the  local  Rleaann  problea.  The 
aaapllng  procedure  la  due  to  van  dcr  Corput  Kef  (S) 
and  generates  equldlatrlbuted  aequencea  Ref  (6).  The 
technique  ellalnatea  nuiierlccl  dlffualon  and  flow 
dlacontlnultlea  ere  computed  with  Infinite  reaolutlon, 
l.e.,  there  la  no  aacarlng. 

The  Godunov  method  la  a  finite  difference  aethod 
which  coaputea  the  Rleaann  problea  aolutlona  aa  a 
flrat  atep,  and  ualng  the  fluxea  thua  obtained, 
advancea  In  tlae  by  aolvlng  the  flrat-order  finite 
difference  approxiutlon  of  the  Euler  equatlona. 

Coding  details  for  the  RCM,  and  the  Godunov 
aethod  are  given  In  Ref a  (4,7)  and  Refs  (7,8) 
reapectlvely.  In  the  following  aeetlon,  ezaaples  of 
the  applleatlone  of  theae  technlquea  to  wave  rotor 
devlcca  are  given. 

EXAMPLE  APPLICATIONS 

Fig  1  ahowa  a  wave  dlagraa  (a  one-dlacnalonal 
apace-tlae  plot)  for  a  wave  rotor  device  configured  to 
operate  aa  a  turbine  alone.  Air  at  a  preaaure  higher 
than  aablent  entera  the  rotor  through  the  Inlet  port, 
generating  a  ahock  wave  with  a  alower  moving  contact 
dlacontlnulty  behind  It  at  point  'a*.  The  ahock 
reflects  off  a  wall  boundary  at  point  'b',  crosaea  the 
Incoming  Interface  at  point  'c'  (bringing  It  to  a 
halt),  and  reaches  the  Inlet  aide  again  at  point  'd', 
whereupon  the  Inlet  port  la  closed.  The  gases  In  the 
rotor  passages  ere  now  In  an  essentially  quiescent, 
high  pressure  state.  At  point  'a',  the  outlet  port  la 
opened  to  a  lower  preaaure,  a  rarefaction  fan  being 
generated  In  the  process.  The  expansion  waves  travel 
to  the  left,  reflect  off  the  solid  boundary  and 
propagate  back  to  the  outlet  side.  The  port  la  closed 
when  conditions  In  the  passages  are  essentially  the 
sasM  aa  at  the  beginning  of  the  cycle. 


Fig.  1  Wave  Diagram  Computed  by  1-D  Random  Choice 
Method. 

S*Shock;  RS*Reflected  Shock;  R-Rarefactlon 
Fan;  RR-Reflected  Rarefaction;  l-Interface 


The  entire  wave  diagram  was  generated  by  the 
one-dlmenslonal  sampling  method  with  properly 
Implemented  boundary  conditions.  Fig  2  shows  a 
sequence  of  the  propagation  of  the  shock  wave  and 
contact  dlacontlnulty  In  time.  The  perfect  reaolutlon 
of  Che  discontinuities  la  noteworthy,  even  when  a  very 
small  density  change  occurs  across  the  Interface,  as 
In  this  case. 


Fig.  2  Sequence  Showing  Shock  and  Interface  Movement. 

S  -  Shock;  I  -  Interface;  RS  -  Reflected 

Shock;  N  -  Times cep 

The  transient  process  of  the  rotor  passage  end 
opening  or  closing  creates  losses  which  affect  the 
performance  of  the  device.  The  gradual  (as  opposed  to 
Instantaneous)  opening/closing  of  the  passages  Is 
essentially  a  Cwo-dlmenalonal  flow  phenomenon  and  Is 
modelled  using  the  2-D  Godunov  code.  Fig  3  shows  a 
sequence  of  pressure  contours  for  Che  gradual  opening 
case  which  clearly  ahowa  the  significant  effect  the 
dynamics  of  the  passage  opening  has  on  the  flow 
pattern  In  wave  rotor  devices.  The  actual  ahock 
formation  for  the  modelled  case  Is  seen  to  occur 
approximately  halfway  Into  the  passage  as  opposed  to 
the  Instantaneous  formation  generally  assumed  for  the 
Ideal  case. 

DISCUSSION 

The  two  techniques  described  here  are  powerful 
tools  for  the  design  of  wave  rotor  devices.  The 
example  application  of  the  one-dlmenslonal  Random 
Choice  method  deals  with  a  fairly  elementary  wave 
dlagraa,  but  serves  to  Illustrate  Its  potential  and 
suitability  for  designing  complex  gas  turbine  engine 
cycle  type  wave  diagrams.  The  Godunov  method  la  more 
amenable  to  extension  to  multidimensional  analysis 
required  for  loss  calculations,  while  fulfilling  the 
requirements  outlined  In  the  Introduction. 
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Pig.  3  Contour  Plots  Generated  by  2-D  Godunov  Method 
for  Gradual  Opening  of  Have  Rotor  Passage  End. 
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Local  Cell  Orientation  Method. 


Illustration  of  the  method  for  the  solution  with  oblique  shock 


INTRODUCTION 


The  Finite  Volume  Schemes  for  numerical  solution  of  the  Euler  equa¬ 
tions  can  be  split  into  two  main  steps: 

a) Calculation  of  the  fluxes  on  the  cell  edges; 

b)  Updating  the  values' at  the  center  of  the  cell  ,  using  the 

calculated  fluxes,  to  satisfy  the  Euler  conservation  laws. 

Usually  the  accuracy  of  the  flux  calculation  determine  the  accuracy 
of  the  numerical  modeling  Ref.1.  In  the  region  of  the  shock  wave  the 
accuracy  of  the  flux  calculation  is  significantly  reduced  and  in  some 
methods  special  care  should  be  taken  in  order  to  reduce  the  pre-shock 
and  after-shock  oscillations .  When  the  upwind  methods  are  used  the 
fluxes  are  usually  calculated  using  the  one  dimensional  wave  propagation 
problem  (i.e.  Riemann  problem).  In  this  case  the  assumption  of  the 
one-dimensionality  lead's  to  large  computational  errors  and  significant 
shock  smearing  in  the  region  of  the  oblique  shock. 

In  the  current  study  it  will  be  demonstrated  that  the  accuracy  of 
the  numerical  modeling  could  be  significantly  improved  by  proper  Local 
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Cell  Orientation  (SELCO  Method).  We  applyed  the  SELCO  Method  only 
with  the  Godunov  scheme  for  the  solution  of  the  Euler  equations,  but 
the  idea  of  the  method  is  general  enough  to  be  used  with  other  finite 
volume  methods. 


THE  SELCO  METHOD. 


For  the  Godunov  method  the  numerical  fluxes  are  calculated  using 
the  solution  of  the  Riemann  problem.  First  the  solution  is  assumed  to  be 
piece-wise  constant  in  every  cell.  Then,  to  calculate  the  flux  for  every 
cell  edge  the  Riemann  problem  is  solved  between  the  cells  adjacent  to 
the  edge.  Solution  of  the  Riemann  problem  is  obtained  in  the  direction 
normal  to  the  cell's  edge.  The  equations  are  then  integrated  in  the 
each  cell  for  one  time  step  using  the  finite  difference  approximation  of 
the  Euler  equations  and  the  fluxes  at  the  cell’s  edges.  A  more  complete 
description  of  the  Godunov  method  can  be  found  in  Ref.1. 

In  the  two  dimensional  Godunov  method  it  is  assumed  that  the  fluxes 
could  be  determined  using  the  solution  of  the  one  dimensional  Riemann 
problem  if  the  CFL  number  used  is  smaller  than  0.5.  In  the  following 
example  it  will  be  demonstrated  that  that  assumption  lead's  to  consider¬ 
able  errors  for  the  solutions  with  oblique  shock  waves. 

As  a  test  case  the  supersonic  flow  over  the  22.3“  wedge  was  simu¬ 
lated.  The  grid  and  the  computational  domain  is  shown  on  the  Figure  1. 
The  system  of  the  two  dimensional  Euler  equations  as  in  Ref.1,  was 
solved  by  time  marshing  in  the  computational  domain  shown  in  Fig.l  for 
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time  t-^oo  with  the  following  conditions  at  the  domain  boundaries: 

a)  Inflow  of  the  air  at  M=2.2  along  segment  1-2; 

b)  Outflow  along  segment  3-4; 

(the  flow  along  this  segment  will  be  supersonic  and 
the  continuation  condition  is  applied  to  all  parameters 
at  the  downstream  boundary) 


c)  Solid  wall  along  segments  1-3  and  2-4. 


In  the  Figure  2  (a,b,c,d)  results  of  the  flow  simulation  using  the 
Godunov  method  in  the  computational  domain  shown  on  Fig.l  (squares) 
are  compared  with  the  analytical  solution  for  the  supersonic  flow  over 
the  wedge  (solid  lines)  Ref. 2.  Comparison  is  done  for  the  data  on  the 
lower  surface  of  the  channel  shown  on  the  Fig.l.  It  can  be  concluded 
from  this  comparison  that  only  pressure  coefficient  is  calculated  accu- 
ratly  by  the  Godunov  method.  Especially  notable  the  error  in  entropy. 
Exact  entropy  change  on  the  shock  wave  is  2  times  smaller  than  that 
predicted  by  the  Godunov  method.  The  isomach  lines  are  shown  for  the 
same  simulated  case,  in  Figure  2  (e).  The  shock  is  smeared  over  the 
large  area  of  the  mesh  (up  to  5-6  points),  but  the  shock’s  angle  is 
simulated  correctly . 

In  order  to  study  the  source  of  errors  for  this  type  of  problems  the 
supersonic  flow  over  the  wedges  with  different  angles  was  simulated.  It 
was  found  that  the  pressure  was  predicted  accuratly  for  the  all  simu¬ 
lated  cases.  At  the  same  time  the  error  in  entropy  increases  when  the 
shock  wave  obliquity  increases.  Dependence  of  the  accuracy  of  the 
shock  simulation  on  the  shock  obliquity  is  not  surprising  since  the  basic 
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assumption  about  independance  of  the  flux  calculation  for  the  intersect¬ 
ing  edges  of  the  cell,  for  t  <  0.5*CFL,  became  incorrect  in  the  vicinity 
of  the  shock.  Spliting  the  flux  calculation  for  each  cell  on  four  one 
dimensional  Riemann  problems  for  each  of  the  cell's  edges  would  not 
lead  to  the  large  errors  in  the  regions  of  monotonic  flow.  It  will  also 
give  accurate  approximation  in  the  region  of  the  oblique  shock  if  the 
shock  wave  is  parallel  to  one  of  the  cell  edges.  The  last  statement  will 
be  illustrated  by  the  following  numerical  example. 

Let's  solve  the  problem  of  the  supersonic  flow  with  M  -2.2  over  the 
22.3°  wedge,  as  in  previous  example,  on  the  grid  shown  in  Figure  3. 
The  grid  in  Fig.l  is  skewed  uniformly  so  that  vertical  lines  of  the  grid 
are  forming  a  51.7°  angle  with  the  positive  direction  of  the  x-axis.  An 
angle  of  51.7°  corresponds  to  the  shock  wave  angle  calculated  analyt¬ 
ically  Ref. 2.  Because  the  grid  is  skewed  the  oblique  shock  waves  will 
be  parallel  to  the  vertical  grid  lines  of  the  mesh  shown  in  Fig. 3. 
Results  of  this  test  case,  in  the  form  of  the  distribution  of  pressure 
coefficient,  entropy,  density  and  velocity,  are  compared  to  the  analyt¬ 
ically  calculated  values  (solid  lines)  in  Figure  4  (a,b,c,d).  In  Fig. 4 
(e),  the  isomach  lines  for  this  simulation  are  shown.  It  can  be  con¬ 
cluded  from  the  results  presented  in  Fig. 4,  that  the  flow  is  modeled  is 
this  case  with  very  high  accuracy.  Shock  wave  is  resolved  on  one  grid 
cell  and  entropy  jump  is  calculated  exactly.  All  that  is  result  of  the 
proper  cell  orientation  towards  the  shock  wave. 

Skewing  the  entire  grid  can  help  to  accuratly  resolve  only  one 
shock  wave,  and  only  on  the  condition  that  the  shock  wave  angle  is 
known  prior  to  the  solution,  so  it  could  not  be  applied  effectivly.  But 
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if  skewing  of  the  cell  can  be  done  locally  in  the  region  of  the  oblique 
shock  it  can  improve  the  accuracy  of  the  shock  wave  simulation. 

The  proposed  SELCO  method  does  the  cell  reorientation  locally.  The 
method  consist  from  the  following  steps; 

1 )  Integration  of  the  Euler  equations  by  the  Godunov  method; 

2) Define  the  approximate  shock  location  and  the  shock  angle  using 

the  expressions  for  oblique  shocks  Ref. 2. 

(Based  on  the  fact  that  pressure  is  calculated  accuratly 
by  the  Godunov  method); 

3)  Rotate  the  cell  edges  that  are  directed  along  the  shock  about 
their  middle  points  on  the  shock  wave  angle.  So,  after  rotation, 
two  of  the  cell  edges  will  be  parallel  locally  to  the  shock  wave 
surface  (see  Figure  5); 

4) Calculate  the  fluxes  on  the  new  edges  of  the  cell; 

5)  Integrate  the  Euler  equations  in  the  new  cell. 

All  these  additional  steps  do  not  add  much  of  computational  work, 
because  the  cell  reorientation  should  be  done  only  in  the  vicinity  of  the 
shock  wave.  If  the  shock  could  be  resolved  on  one  grid  point,  only  one 
cell  would  be  transformed.  Our  experience  has  shown  that  there  is  no 
need  to  extend  the  SELCO  procedure  on  the  cells  ahead  and  behind  the 
shock  wave  surface. 

In  Figure  5  (a,b,c,d)  results  are  shown  for  the  same  flow  condition 
as  in  previous  cases:  supersonic  flow  with  M  =2.2  over  the  wedge  of 
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22.3°.  Calculations  were  done  on  the  grid  shown  in  Figure  1  and  in 
this  case  the  SELCO  method  was  applyed.  It  can  be  concluded  from  the 
results  presented  in  Fig.  5  that  the  accuracy  of  the  shock  wave  model¬ 
ing  using  the  SELCO  method  approach  that  demonstrated  in  Figure  4  for 
the  completely  skewed  grid,  and  is  superior  to  the  accuracy  of  the 
standard  Godunov  method.  The  isomach  lines  for  the  case  where  the 
SELCO  method  was  applied  are  presented  in  Figure  5  (e).  The  shock 
tickness  in  this  graph  is  minimal  and  much  more  improved  compared  to 
the  simulation  by  the  original  Godunov  method  shown  in  Figure  2  (e). 


CONCLUSIONS. 


It  has  been  demonstrated  that  inaccurate  modeling  of  the  oblique 
shock  waves  produced  by  the  Godunov  method  is  the  result  of  the 
obliqueness  of  the  shock  wave  with  respect  to  the  edges  of  the  cells  of 
the  computational  grid  covering  domain  of  integration.  It  is  also  shown 
that  only  when  the  shock  surface  is  parallel  to  the  two  opposite  edges 
of  the  cell  the  oblique  shock  can  be  accuratly  calculated. 

A  new  method  of  the  Local  Cell  Orientation  {SELCO  Method) 
is  proposed  in  order  to  allow  local  reorientation  of  the  cells  in  the 
vicinity  of  the  shock  waves.  The  efficiency  of  the  SELCO  method  is 
demonstrated  for  the  simulation  of  the  oblique  shock  waves  in  the 
supersonic  flow  using  Euler  equations. 


Although  the  new  SELCO  method  was  demonstrated  with  the  Godunov 
scheme  it  will  be  effective  in  applications  to  the  other  upwind  methods 
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that  use  the  finite  volume  formulation. 
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FIGURE  CAPTIONS. 

1.  The  Computational  Domain. 

2.  Computation  by  the  Godunov  method.M^=2.2. 

Wedge  Angle  is  22.3°.  Continuous  lines 
show  the  exact  solution  for  -  a,b,c  and  d. 

a)  Surface  Pressure  Coefficient. 

b)  Surface  Entropy. 

c)  Surface  Density. 

d)  Surface  Mach  Number. 

e)  Isomach  Lines. 

3.  The  Computational  Domain  with  the  Skewed  Grid. 

4.  Computation  by  the  Godunov  Method  on  the  Skewed  Grid  Shown 
in  Figure  3.  M.^=2.2.  Wedge  Angle  is  22.3“. 

Continuous  lines  show  the  exact  solution  for  -  a,b,c  and  d. 

a)  Surface  Pressure  Coefficient. 

b)  Surface  Entropy. 

c)  Surface  Density. 

d)  Surface  Mach  Number. 

e)  Isomach  Lines. 

5.  The  Local  Cell  Reorientation  by  the  SELCO  Method. 

6.  Computation  by  the  SELCO  Method .M_^=2. 2. 

Wedge  Angle  is  22.3“.  Continuous  lines 
show  the  exact  solution  for  -  a,b,c  and  d. 

a)  Surface  Pressure  Coefficient. 

b)  Surface  Entropy. 

c)  Surface  Density. 

d)  Surface  Mach  Number. 

e)  Isomach  Lines. 
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a)  Surface  Pressure  Coefficient.  c)  Surface  Density. 

4.  Computation  by  the  Godunov  Method  on  the  Skewed  Grid  Shown 


In  Figure  3.  M.^=2.2.  Wedge  Angle  is  22.3*. 

Continuous  lines  show  the  exact  solution  for  -  a,b,c  and  d, 


Wedge  Angle  is  22.3*.  Continuous  lines 
show  the  exact  solution  for  -  a,b,c  and  d. 


